/***********************************************************
*  --- OpenSURF ---                                        *
*  This library is distributed under the GNU GPL. Please   *
*  contact chris.evans@irisys.co.uk for more information.  *
*                                                          *
*  C. Evans, Research Into Robust Visual Features,         *
*  MSc University of Bristol, 2008.                        *
*                                                          *
************************************************************/

#include "utils.h"

#include "integral.h"

using namespace std;

//! Computes the integral image of image img.  Assumes source image to be a
//! 32-bit floating point.  Returns IplImage of 32-bit float form.
IplImage *Integral(IplImage *source)
{
    // convert the image to single channel 32f
    IplImage *img = getGray(source);
    IplImage *int_img = cvCreateImage(cvGetSize(img), IPL_DEPTH_32F, 1);

    // set up variables for data access
    int height = img->height;
    int width = img->width;
    int step = img->widthStep/sizeof(float);
    float *data   = (float *) img->imageData;
    float *i_data = (float *) int_img->imageData;

    // first row only
    float rs = 0.0f;
    for (int j=0; j<width; j++)
    {
        rs += data[j];
        i_data[j] = rs;
    }

    // remaining cells are sum above and to the left
    for (int i=1; i<height; ++i)
    {
        rs = 0.0f;
        for (int j=0; j<width; ++j)
        {
            rs += data[i*step+j];
            i_data[i*step+j] = rs + i_data[(i-1)*step+j];
        }
    }

    // release the gray image
    cvReleaseImage(&img);

    // return the integral image
    return int_img;
}

//-------------------------------------------------------

//! Computes the sum of pixels within the rectangle specified by the top-left start
//! co-ordinate and size
inline float BoxIntegral(IplImage *img, int row, int col, int rows, int cols)
{
    float *data = (float *) img->imageData;
    int step = img->widthStep/sizeof(float);

    // The subtraction by one for row/col is because row/col is inclusive.
    int r1 = std::min(row,          img->height) - 1;
    int c1 = std::min(col,          img->width)  - 1;
    int r2 = std::min(row + rows,   img->height) - 1;
    int c2 = std::min(col + cols,   img->width)  - 1;

    float A(0.0f), B(0.0f), C(0.0f), D(0.0f);
    if (r1 >= 0 && c1 >= 0) A = data[r1 * step + c1];
    if (r1 >= 0 && c2 >= 0) B = data[r1 * step + c2];
    if (r2 >= 0 && c1 >= 0) C = data[r2 * step + c1];
    if (r2 >= 0 && c2 >= 0) D = data[r2 * step + c2];

    return std::max(0.f, A - B - C + D);
}